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Abstract 


Possibly the best set of data for earthquake excitation of soils exists for the test site operated by the 
Electric Power Research Institute (EPRI) and the Taiwan Power Company at Lotung Taiwan. At 
this site, two locations are instrumented with three-component accelerometers at depths of 47, 17, 
11, 6 meters, and at the surface. One array is in the free-field while the other is adjacent to a one- 
quarter scale nuclear containment vessel. The site is also well instrumented with piezometers at 
various depths and locations. The simplified soil profile consists of 30 to 35 m of silty sand and 
sandy silt with some gravel, overlaying a thick clay and silt deposit. The water table is within half 
a meter of the ground surface. This area is seismically active, and strong shaking generated by 


many earthquakes exhibiting a wide range of magnitudes have been recorded since 1986. 


This report summarizes the data and signal processing that was done to the EPRI Lotung data at 
the Colorado School of Mines. The over 2000 files were organized into MATLAB experiments 
by event. The provided acceleration data were carefully double integrated to yield velocity and 
displacement time history records. The data are now in proper format to begin system identifica- 


tion analysis. 
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CHAPTER 1 - INTRODUCTION 


1.1 Background 


There are many unanswered questions of interest to the geotechnical community concerning the 
behavior of soils subjected to earthquake excitation. Foremost among them are questions 
concerning the strain-dependent non-linear behavior of soils, and soil-structure interaction. In an 
attempt to gain further rational understanding of these problems the Electric Power Research 
Institute (EPRI), and Colorado School of Mines (CSM), and the Structures Division of the National 
Institute of Standards and Technology (NIST) formed a cooperative research team to evaluate 


ground motion time histories recorded at the Lotung site, Taiwan. 


Much of the necessary data is to be made available by EPRI from the Lotung site. The Lotung 
strong motion data set are an extremely unique set of data. The completeness of this input-output 
data set makes it ideal for analysis using system identification (SI) methods. Data from the 
Wildlife Site, Imperial Valley CA, and Treasure Island are freely available. This report lays out 
the work undertaken by the P.I. at CSM for 1994, with funding provided by NIST under Award 
Number 60NANB4D1677. 


1.2 Why Use System Identification? 


An important goal for earthquake engineering is the ability to estimate soil properties without 
intruding into the soil mass. For the engineer interested in seismic behavior of soils, the dynamic 
properties of the soil are of interest, particularly large-strain properties. The archetypal large strain 
field excitation is earthquake strong motion. Ideally, both ground motions into the soil layer of 
interest and on the surface above the layer are recorded, as illustrated by the cartoon in Fig. 1.1. 
Given this known input propagating upward from depth, and the output at the top of the soil 
column, the behavior of the soil can be modelled by inverse theory. If a suitable model is chosen 
to represent the system of interest, the estimated model parameters will correspond to important 
mechanics parameters of the system, such as damping, natural frequency, and stiffness. This 


estimation of parameters is commonly known as system identification (SI). 
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Fig. 1.1 - Configuration of the System Identification Method. 


The traditional method of geotechnical analysis of dynamic soil motions is through the Fourier 
transform. However, serious problems arise when this method is applied to short data streams, and 
to signals changing through time — non-stationary signals. This study was undertaken to show the 
effectiveness of a different type of model, a parametric model commonly used in automatic control 
and geophysics, which avoids many of the limitations inherent with calculating the system transfer 
function by Fourier techniques. An important aspect of particular parametric models is the 
theoretical link between the estimated system parameters and the mechanical parameters of a 
lumped-mass oscillator. The parametric model allows estimates of system dynamic properties to 


be to be made if an input-output data set is available. 


1.3 Scope 


The purpose of this report is to set the stage for future detailed SI analyses of the Lotung site. To 
this end the Lotung site itself will be introduced through geological, seismological, and 
geotechnical description. The method of data processing used to prepare the ground acceleration 
and pore water pressure time records will be discussed in detail. Finally a brief example of SI 


analysis on an individual record will be presented to illustrate the method to be used. 


CHAPTER 2 - THE LOTUNG SITE, TAIWAN 


2.1 Introduction and Geography 


With the growth of the use of nuclear-powered generating plants in the 1970’s, many safety related 
questions about the seismic performance of these plants arose. In the early 1980’s, EPRI and the 
Taiwan Power Co. (TPC) constructed two scale models (1/4 and 1/12 scale) of a nuclear 
containment structure near Lotung, Taiwan. This is a very seismically active area in northeast 
Taiwan (see Fig. 2.1). The site and structures were elaborately instrumented so that soil and 
structural response, and soil-structure interaction, to earthquakes could be carefully studied (Tang 


et al., 1989; Liu and Yeh, 1985). 
ee Instrumentation 


The soil instrumentation includes a three-arm surface array, as shown in Fig. 2.2. The arms radiate 
approximately 47 m from the 1/4 scale containment structure. In addition, there are two downhole 
arrays of accelerometers extending to a depth of 47 m, as shown in Fig. 2.2. The surface 
accelerometers are triaxial force-balance units (Kinemetrics FBA-13) oriented in the N-S, E-W, 
-and vertical directions. The downhole arrays (DHA and DHB) are modified Kinemetrics FBA- 
13H units oriented in the N-S, E-W, and vertical directions. DHA is located 3 m from the 
containment vessel and DHB is located 47 m from the structure, allowing identification of the 
effects of the structure on soil response. The downhole instruments are located at depths of 6, 11, 
17, and 47 m. The simplified soil profile consists of 30-35 m of silty sand and sandy silt with some 


gravel, above clayey silt and silty clay. The water table is within half a meter of the ground surface. 


PACIFIC 
OCEAN - 


SMART 1 
ARRAY 


Fig. 2.1 Location of Lotung Large-Scale Seismic Test Site (LSST). EPRI (1989) 
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gS) Site Characterization 


The basic geology of the LSST site is summarized by Wen and Yeh (1984) and Tang (1987). The 
area consists of a recent alluvium layer 40 to 50 m thick overlying a Pleistocene formation that 
varies from 150 to 500 m in thickness. Underlying the Pleistocene material is a Miocene basement 
rock. Characteristic geological profiles also showing compressional wave velocities are shown in 
Fig. 2.3. Example soil profiles are shown on Fig. 2.4. The locations of the boreholes from which 


the profiles were constructed are shown on Fig. 2.5. 


Five stages of laboratory testing programs were performed at the LSST during specific phases of 
the project to determine engineering properties of the soil. A summary of the tests performed is 


included here and the specific references for the test results are given below. 


1. 1984 Jong Shing Boring Services (JSBS) Laboratory Testing Program 


Index properties: 
Soil classification 
Grain size analyses 
Moisture contents 
Specific gravity 
Dry density and void ratio 
Atterberg limits 
Engineering properties: 
Direct shear tests 
Triaxial shear tests 
Unconfined compression tests 


2. 1987 National Taiwan University (NTU) Laboratory Testing Program 


Index properties: 
Grain size analyses 
Moisture content 
Specific gravity 
Dry density 
Atterberg limits 
Engineering properties: 
Uniaxial - load/unload and cyclic loading tests 
Triaxial - compression, extension and cyclic loading tests 
Resonant column 
Hydrostatic - load/unload and cyclic loading tests 
Compaction tests 


3. 1987 University of California at Davis (UCD) Laboratory Testing Program 


Engineering properties: 
Triaxial shear loading 
Cyclic triaxial liquefaction testing 


This phase of the laboratory testing was performed in conjunction with the installation of pore 
water pressure transducers at the site. Pore water pressures were monitored in the samples during 


these laboratory tests. 


4. 1989 NTU Laboratory Testing Program 


Engineering properties: 
One dimensional rebound 
Resonant column 
Cyclic triaxial liquefaction 
Permeability 
Cyclic triaxial modulus 
These additional tests were performed by NTU for the specific purpose of investigating the form 
of the shear modulus versus shearing strain and material damping versus shearing strain curves for 


undisturbed soil samples. 


5. 1990 UCD Laboratory Testing Program 
Engineering properties: 

Cyclic triaxial modulus 

Cyclic triaxial liquefaction 

Cyclic simple shear 

One dimensional rebound 

Permeability 
These tests by UCD were independent of the 1989 NTU tests and were performed to investigate 
discrepancies in modulus and damping data found from analysis of 1987 results for the LSST site. 
In addition, it provided additional data on the cyclic strength and liquefaction properties of the soil. 
Blowcount results from the SPT are shown on Fig. 2.6 for 2 of the 3 drilling and sampling 
programs. The appearance of an occasional layer requiring an excess of 50 blows per foot indicates 
a gravelly soil zone that occurs in discrete lenses rather than as a consistent layer. Results of CPT 


soundings are shown in Fig. 2.7 and 2.8. Occasional spikes in tip resistance values confirm the 


presence of the gravelly lenses. 
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After Tang (1) 
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Fig. 2.3 Geological profiles for the Lotung site. EPRI (1989) 
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Fig. 2.4 Detailed soil profiles developed during the WECC field program. EPRI (1989) 
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Fig. 2.5 Location of boreholes at the Lotung site, and location of soil profiles. 
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Fig. 2.6 SPT blowcounts recorded at Lotung test site. 


Anderson (1993) 


25 


SILTY SAND 


CLAYEY SILT 
SAND OR SILT 


CLAYEY SILT 
SILTY SAND 


SILTY CLAY 


SILTY SAND 


SILTY CLAY 
SAND 


SILTY SAND 


SILTY CLAY 


SILTY SAND 


SAND 


0 50 100 =©180 += 200-~Ss«0.00 0.05 0.10 
TIP RESISTANCE ( ksc ) FRICTION RATIO SOIL PROFILE 


Fig. 2.7 CPT tip resistance, friction ratio and estimated soil profile at test hole #1. | Anderson (1993) 
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Fig. 2.8 CPT tip resistance, friction ratio and estimated soil profile at test hole #2. 


2.4 Seismology 


A summary of the properties of measured temblors is given in Table 2.1. The epicenters for these 
events is shown on the map presented in Fig. 2.9. For the initial interpretation work it was decided 
to concentrate on the events for which accurate pore pressure records exist and show some increase 
in pore pressure during strong shaking. Events 12, 16, and 17 meet this criteria. Events 12 and 16 
were major events and have been discussed in detail (e.g. EPRI, 1989; Chang et al., 1991; 
Anderson, 1993). Event 17 is an aftershock of event 16 which generated enough pore pressure to 
register on the recording equipment. The time histories for acceleration, velocity, and 


displacement for Event 16 are shown in Fig. 3.1, with a full velocity records given in Appendix A. 


2.5 Pore Pressure Generation 


Over the history of the Lotung experiment, Over 30 pore water pressure transducers were installed 
at the site (Shen et al., 1989). The location of these sensors in relation to the three-arm surface 
array is Shown in Fig. 2.10. The soil conditions at several pertinent locations are given in Table 
2.2. As reported by Shen et al. (1989), most of the sensors ailed due to mechanical problems. 


However, several remained in operating condition and were triggered by events 12, 16, and 17. 


A typical pore pressure record is shown in Fig. 2.11 in relation to the acceleration time history. 
Information as to in situ pore pressure and increase for each sensor for events 12, 16, and 17 are 
given in Tables 2.3, 2.4, and 2.5, respectively. The time histories from events 12, 16, and 17 are 


given in Appendix B. 


LSST 
LSST 
LSST 
LSST 
LSST 
LSST 
LSST 
LSST 
LSST 9 
LSST 10 
LSST 11 
LSST 12 
LSST 13 
LSST 14 
LSST 15 
LSST 16 
ESSie 4 
LSST 18 
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Event No. 


Peak Acceleration 


Focal 
Depth 


Date Magnitude | Epicentral 
ee [Baan 

9/20/85 
10/26/85 
11/ 7/85 
1/16/86 

3/29/86 

4/ 8/86 

5/20/86 

5/20/86 

7/11/86 

7/16/86 

7/17/86 

7/30/86 

7/30/86 

8/ 5/86 

11/14/86 
11/14/86 
11/14/86 
11/15/86 


Table 2.1 Properties of the recorded LSST series of earthquakes. 
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Fig. 2.9 Locations of the epicenters for LSST seismic events 1 through 18. Tang, et al (1992) 
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Fig. 2.10 Location of the pore water pressure sensors ant the Lotung test site. 
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EPRI (1989) 


PA-1 w(%) PF-1 w(%) PF-5 w(%) 


F = 15.2 7.8 F = 10.9 12.4 F = 92.7 SBiae} 
SM+G SM-SC SM 
CH 1 CHES CH 9 
D = 5.06M D = 3.25M DO = 12.0M 
T= ML T = ML T = ML 
F =68.5% 37.8 F = 56.3% 40 F = 65.5% 31.6 
R=9 R= 18 ReanS 
F = 14.6 F = 50.2 2120 F = 26.9 21.9 
ML SM 

PN1-0 PN1-1 PF-2 
F = 10.9 254 F = 15.2 7.8 F = 92.4 40 
SM+G cL CL-ML 
CH ll CH 12 CH 6 
D = 3.16M D = 6.03M D = 6.05M 
T = ML T= SM T = ML 
F = 92.4% 40 F = 43.8% 37.8 F = 50.2% 68 
R= 5.5 R = 25.5 R=4 
fa= 5022 27.0 F = 14.6 24.3 F = 57.0 BYTE: 
ML-CL CL—ML 

PN2-1 PA-3’ PN3-1 
FS jN@es Ley F = 8.6 Ut F = 47.6 21.4 
SM SM+G Sm-SC 
CH 18 CH 21 CH 24 
D = 6.30M D = 5.10M DO = 6.30M 
T = SM T = SM+G T= SM 
F = 19.5% 14.6 F = 32.0% des: F = 39.9% BYLT/ 
R = 16.5 (23) R = 17 (20.4) R = 12 (12.4) 
F = 42.7 28.3 F = 64.9 lire 
SM x ML 

PN3-2 PN2-2' 
Fa=- 97.5 BOSE F = 10.1 Lye) 
CL SM 
CH 29 W: water content CH 23 
D = 11.0M F: fine-grained fraction D = 8.00M 
T = ML R: bu/Ug (%) T = SM 
F = 68.6% Se, OD: depth of sensor F = 19.5% 14.6 
Re=—> T: soil type R = 4.3 
fea 17, 35.8 Fi = 4227 28.3 
SM SM 


Table 2.2 Soil conditions at the locations of pore water sensors. | EPRI (1989) 
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Fig. 2.11 Typical pore pressure record from event 16, LSST. Depth = 6.38 m; initial hydrostatic 
pressure = 72.6 kPa; increase in pore water pressure = 8.68 kPa. 


Sensor 
No. 


PA-1 
PF-1 
PF-2 
PF-5 
PN1-0 
PN1-1 
PN1-4 
PN2-1 
PA-3' 
PN3-1 
PN3-2 


Notes: 


Channel Depth, h 
No. (m) 
1 5.06 
5 3.25 
6 6.05 
9 12.00 
11 BENG 
12 6.03 
15 5.53 
18 6.30 
21 5.10 
24 6.38 
29 11.00 


u, (2) 
(kPa) 


54.60 
32.80 
55.70 
121.85 
B92 72 
61.30 
62.34 
62.12 
56.15 
72.60 
116.60 


(3) 
(kPa) 


59.90 
38.97 
57.67 
126.14 
42.12 


123.41 


Event LSST12 


Au“? 
(kPa) 


5.30 
6.17 
1.97 
4.29 
2.40 
17.05 
Seg 
11.24 
13.73 
8.68 
6.81 


(1) This table was modified from reference (2). 
(2) u, = Hydrostatic pressure based on initial readings of the 


recordings. 


(3) u, = Recorded peak pore water pressure. 
(4) Au = u, - uy = Peak induced pore water pressure. 


Table 2.3 Peak pore water pressures, event 12. 


Tang, et al (1992) 


File 
Name 


E12C0Ol. 
E12C05. 
E12C06. 
E12C09. 
E12C11 
E12C12. 
E12C15. 
E12C18. 
E12C21. 
E12C24. 
EV2ZC297 


PWP 
PWP 
PWP 
PWP 
. PWP 
PWP 
PWP 
PWP 
PWP 


PWP . 
PWP 


Sensor 
No. 


PF-8 
PN2-1 
PA-3' 
PN2-2' 
PN3-1 


Notes: (1) 
(2) 


(3) 
(4) 


Channel Depth, h Event LSST16 


No. ~ (m) Ae uo? ree File 
(kPa) (kPa) (kPa) Name 

17 15.00 114.60 135.35 20.75 E16C17.PwP 

18 6.30 66.90 83.60 16.70 E16C18.PWP 

21 5.10 56.30 68.74 12.44 E16C21.PWP 

23 8.00 97.20 106.61 9.41 E16C23.PWP 

24 6.38 80.10 90.44 10.34 E16C24.PWP 


This table was modified from reference (2). 

u, = Hydrostatic pressure base on initial readings of the 
recordings. 

u, = Recorded peak pore water pressure. 

4u = u, - u, = Peak induced pore water pressure. 


Table 2.4 Peak pore water pressures, event 16. Tang, et al (1992) 
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- Sensor Channel Depth, h Event LSST17 
No. Bor (m) u, ©? (3) Au’ File 
(kPa) (kPa) (kPa) Name 
PF-8 ee, 15.00 114.60 128.73 14.13 E17C17.PWP 
PN2-1 18 6.30 66.90 72.77 5.87 E17C18 . PWP 
PA-3' 21 5.10 56.30 61.08 4.78 E17C21.PWP 
PN2-2' 23 8.00 97.20 101.92 4.72 E17C23.PWP 
PN3-1 24 6.38 80.10 84.48 4.38 E17C24.PWwP 


Notes: (1) This table was modified from reference (2). 
(2) u, = Hydrostatic pressure based on initial readings of 
event LSST16 recordings. 
(3) u, = Recorded peak pore water pressure. 
(4) Au =u, - u, = Peak induced pore water pressure. 


Table 2.5 Peak pore water pressures, event 17. Tang, et al (1992) 


Chapter 3 - Procedures for Data Processing and System Identification 


3.1 Data Processing 
3.1.1 Cataloging 


The EPRI-supplied data consists of 10 pc-format floppy disks of data. The data supplies ground 
motions, structural motions, dynamic earth pressure, and ancillary data for 18 earthquake events 
(Tang and Tang, 1992). The data set consists of 2,103 individual files. Pore water pressure values 
are available for events 12, 16, and 17 only. As an example of the state of the “raw” supplied data, 
the complete set of “raw” ground motions and pore pressure time histories for event 16 is presented 


in Appendices A and B, respectively. 


The first step of the data processing was to enter the relevant data files into MATLAB and group 
them into corresponding event files. Files event01 through event18 were created. A simplified 
naming convention was developed for each data record. The first letter of the name is an a, v, or 
d, for acceleration, velocity, and displacement, respectively. The next three letters are dha or dhb 
for downhole array a or b. The following number refers to the depth in meters - 0, 6, 11, 17, and 
47. The second number, following the underscore refers to the event number. The last letter - e, 


n, or u, refers to accelerometer orientation - east-west, north-south, or up-down. 


An example of the naming convention is adha47_18u. This record is the acceleration time history 
at a depth of 47 m at downhole array a, event 18. The record gives the vertical ground motion for 


this location and event. 
3.1.2 Filtering, Resampling, and Integration 


The data as received from EPRI is in the form of raw acceleration records, which we processed 
following standard U.S.G.S. method (Converse and Brady, 1992). The acceleration records are 
digitized at a rate of 200 samples per second (ss), for a bandwidth of 100 Hz. This relatively high 
Nyquist frequency causes the event records to be very long, 8,000 data points for each 40 s trace. 
In addition, there is little useful information for our study above 10-15 Hz, and this region would 
be very noisy. Resampling greatly reduces the amount of data to be later analyzed, and eliminates 


time-domain aliasing of the band-passed signals. It was therefore decided to low-pass filter and 
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resample the data at a rate of 25 ss, yielding a Nyquist frequency of 12.5 Hz. Previous work 


indicated that the information of interest would be contained in this band (Glaser, 1993, 1995a.,b). 


The data was resampled using the resample algorithm from the MATLAB Signal Processing 
Toolbox (Krauss et al., 1994). The data is first low-pass filtered using a Kaiser-windowed linear- 
phase FIR filter using ten terms on either side of the timestep in the calculation. The low-pass 
filter is applied both forward and reverse to eliminate phase shift, and the data resampled at 25 Hz. 
The same process was carried out for every strong motion record as well as for the pore water 


pressure records. 


The acceleration records were then processed and integrated to provide velocity records, and the 
velocity records similarly reprocessed and integrated to yield displacement time histories. The 
procedure used to accomplish these transformations is presented in Appendix C. The acceleration 
record is first high-pass filtered at 0.15 Hz using a 4th order bi-directional Butterworth filter to 
remove dc offset and low frequency drift. A best-fit straight line from the arrival of the strong 
motion is then subtracted from the data, followed by the sample mean. The pre-arrival data is then 
set to zero. The processed acceleration data is now integrated using the trapezoidal method 
(Converse and Brady, 1992). The same steps are carried out on this newly formed velocity time 
history to yield the strong motion displacement. A typical set of processed ground motion records 


is shown in Fig. 3.1. 


3.2 System Identification 
3.2.1 Parametric Modelling 


The goal of system identification is to model a system in a manner that provides needed mechanical 
information about that system. The most common techniques have evolved from electrical and 
mechanical engineering, and involve solving the inverse problem for the system transfer function. 
Each method has limitations; in the words of G. E. P. Box, "All models are incorrect, but some are 


more useful than others." 


24 


Raw Acceleration adha0_16eprelim 
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Fig. 3.1 Typical ground motion record set at Lotung, east-west component of the surface of 
downhole array a, event 16. 


2 


The process of inversion allows the estimation of the system response function (filter) if the input 
and output signals are known. A simple model for characterizing a system is as a parametric 
relationship between system input and output. Such a model, referred to as an autoregressive- 


moving average (ARMA) model, is based on discrete time series analysis: 


nb na 
Vy = OY yh ai he tat -( » bx, _;+ Dy or] (1) 
j=0 k=1 


where y; is the actual output data sequence, x, is the input sequence (assume white noise for simple 


spectral estimation), na and nb are the AR and MA orders, respectively, and the subscript is the 
time step counter. The output is seen as a combination of the input history acted upon by the "b" 
coefficients plus the past outputs acted upon by the "a" coefficients. The input series, involving 
the "b" coefficients, is a causal moving average (MA) process (convolutional). The series 
involving weighted past output values ("a" coefficients) is a noncausal autoregressive (AR) 
process. The lengths of the AR and MA processes (model order) must be explicitly chosen so that 


the model best represents the process. 


Applying the shifting theorem to Eq. 1. yields the Fourier transform (Bracewell, 1978) 
Al 1+ a,e" + ase i + sy = Xu( bo + be” + bye? + ay (2) 


where i is J—1 and @ is circular frequency. Applying the Z-transform (Bracewell, 1978), where 


k ki@ : : P : 
z=e ,and rearranging, yields the frequency domain transfer function H,, 


1 2 
be Debi +ibS 7a 
H = ee er eae : i (3) 
@ 12a A520 oer 
The ARMA model is very powerful in that it can easily model sharp drops, sharp peaks, and 
smooth spectral behavior. It is also the most parsimonious estimator (Robinson, 1982), describing 


a complex process with very few parameters calculated from a small length of data. Parametric 


modelling avoids many of the difficulties inherent in the traditional Fourier methods, discussions 
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of which can be found in many books and journals (e.g., Glaser, 1993; Johansson, 1993; Pandit, 
1991). Extensions of this model, e.g., ARMAX, ARX, Box-Jenkins, allow input, system, and 
output noise to be expressly modelled (Ljung, 1987). In particular, the ARX model includes the 


effect of uncertainties and noise as a white noise term. 


The ARMA model has special significance since it can be derived directly from the differential 
equation of motion for an N-degree-of-freedom (DOF) system, with the damping ratio and 
resonant frequency as the model parameters (e.g., Gersch and Luo, 1970). A 2n-2n ARMA model 
is therefore a valid model for a layered soil system, or soil-structure interaction problem. The 
damping ratio and resonant frequency of the N-degree-of-freedom oscillators are contained in the 
2n AR parameters. Phase relations are preserved in the MA parameters. The modal frequencies j, 
percent of critical damping ;, (Ghanem et al., 1991) and power participation factor p; (Pandit, 
1991; Safak, 1988) are calculated from the system poles and residues found from partial-fraction 


expansion of Eq. 2. The modal parameters are defined as 


JA; + 8 | (4) 


OF = — Kt 

E S (5) 
; J; + oe 

P; = oe r,conj (Z;) — Z,conj (7;) (6) 


where 1 = Arg (z,),d = —(0.5) In|z af , Zj is the pole for mode j, 7; is the residue for mode j, and 


At is the digitization rate. 
3.2.2 Adaptive (recursive) model estimation 


_ Traditional methods of system estimation, both parametric and non-parametric, are strictly valid 
only for stationary data. A stationary signal is one whose statistics do not change with time. The 
commonly invoked, loose definition of stationarity requires that the variance of the signal be 


constant over any and all time windows. Inherent in the Fourier transformation of a time series to 
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the frequency domain is the averaging of the signal components over the sampling period T. A 
piece of time is frozen over this period and the assumption made that all time before and after is 
the same, i.e. repeated forever. The energies present at each component frequency are integrated 


over the entire time period T. 


The difficulty with non-stationary signals is that these energies are changing during this period. If 
the frequencies present are changing over this time window, the resulting estimation, regardless of 
method used, will be a smeared average as if all the frequencies with a given energy were active 
throughout the entire period. For weakly non-stationary processes, the effect over a small time 
period is unimportant. If needed, the signal can be cut into relatively stationary sections and 


spectra found using methods specially designed for short data segments, i.e. Burg's method. 


The field of adaptive filtering was formed to model non-stationary processes. As the statistics of 
the signal change through time, the filter "adapts" to the changing variance with new parameters 
that reflect the structure of the system at that point. The predicted value for the next time step can 
be compared with the actual value, and the difference (referred to as innovations) 


Pure 
By, 


where y, is actual output at time rt, and ), is the prediction of output at time ¢ made at time t-J, will 


give a measure of how well the filter is doing its job. The term "innovations" is used because this 


information is new information that can not be predicted by the model at this particular step. 


Autoregressive parameters can be sequentially estimated so that the parameters are adaptive to the 
changing nature of the process (Marple and Lawrence, 1987). The AR parameters are updated 
after each data point, tracking slowly non-stationary signals. A forgetting factor, commonly a 
damped negative exponential, is used so that older data carries less and less weight. A frequency 
domain estimation can be made at any time step by evaluating the AR parameters around the unit 


circle, giving the spectral representation of the behavior of the process at that time. 
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The most popular direct adaptive filter, or process model, is the so-called Kalman filter (Kalman, 
1960; Kalman and Bucy, 1961). Sorenson (1970) points out that the Kalman approach is a direct 
descendant of Gauss's least squares, except now neither the signal nor the noise model must be 
stationary — the state may change from sample point to sample point. Nau and Oliver (1979) state 
that the Kalman filter is based on a dynamic AR model defined by "two concurrent random 


equations of motion": 


vf 
x, — Jat poi + a. (7) 


the AR(p) equations of motion, and the "motions" of the parameters, 


® = ®_,+5, (8) 
where p = number of prior observations utilized, 

He = vector of p prior data observations x Ae perth 

®, = vector of p AR parameters, 

at = Gaussian white noise with 0 mean and variance o7 

by = Gaussian white noise with 0 mean and covariance matrix Q. 


Equation 9 estimates a value of ©; comprised of p previous parameters, through a random walk 
| equation. The estimate uses the weighted p previous data points, and yields a new observation x4 


when added to a new noise value. The least squares solution solves the equations so that the 
innovations (Eq. 7) — new, dynamic information that cannot be predicted — are minimized in a 


least squares sense each time step. 


The theory behind the Kalman filter can be manipulated to yield the system parameters for the case 
where there is no a priori information about the noise, and even when there is no information about 
the input signal. The so-called extended Kalman filter has been very successfully applied to non- 
stationary (and non-linear) estimation problems (Ljung, 1979; Astrom and Eykhoff, 1971). The 
manner of application is actually straight-forward. The Kalman model is constantly updating its 
estimation of the dynamic process by examining the innovations. The dynamics can be due to a 


changing input or noise process, or it can be due to the system itself changing. The effect is a 
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linearization between single time steps, but if the system is changing slowly compared to the time 


step used, the linearization is "invisible" and the non-linear behavior is well modelled. 
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Chapter 4 Parametric Modeling of Lotung Data 


Previous work sponsored by NIST (Glaser, 1995a) demonstrated that the MATLAB software 
package (Mathworks, 1993) is ideally suited for manipulation, processing, and presentation of 
earthquake data. MATLAB is a matrix-based system which evolved from the LINPAC and 
EISPAC libraries commonly used for mainframe FORTRAN numerical analysis. Complex 


numerical problems can be speedily solved without programming in the traditional sense. 


For this study, the standard routines contained in the MATLAB System Identification Toolbox 
(Ljung, 1993, 1987) allowed SI to be used as a tool accessible to the geotechnical engineer. 
Virtually every approach and algorithm encountered in the literature by the author could be 
duplicated rapidly and accurately. When run on an SGI Indy workstation, all aspects of the 


analysis were quick enough to be interactive. 


Analysis begins by determining if the event can be modelled as stationary segments. A recursive 
segmentation scheme, which attempted to break the data into segments with a chosen maximum 
variance (Ljung, 1987), is used. However, it is not possible to determine the “correct” variance a 
priori. In these cases a more direct method is used — the output simulated by the calculated 


system has to accurately model the actual measured output. 


When possible the input-output data record is broken into segments based on a mechanistic 
understanding of the seismic event and soil behavior gained from study of the pore pressure 
behavior. Initially it is assumed that the various segments are basically stationary. If the stationary 
model can not accurately and parsimoniously simulate the segment output, a non-stationary 
recursive model is used. In addition, the appropriateness of the model is checked by insuring 99 
percent confidence in both the whiteness of the residual autocorrelation function and the cross- 


correlation function between the input and output residuals (Bohlin, 1987). 


The stationary SI algorithm uses a least squares estimation for the ARX model. It is necessary to 
estimate the number of parameters to be calculated, which is essentially estimating the degree of 


freedom of the soil system. There is no obvious answer to the degrees of freedom of the system, 
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so several verification techniques are employed to insure that a proper model order is estimated. 
Model order is increased in 2n-2n steps corresponding to an additional degree of freedom each. 
The simulated output of the model is then compared with the actual output for congruence, and the 
fewest parameters needed to accurately characterize the system was chosen as the system model 
order. Examination of the pole and zero plot insures that excessive, overlapping parameters are 
not included (Ljung, 1987). If the segments proved non-stationary, they are analysed using a 


recursive Kalman filter technique that expressly accounts for non-stationarity (Ljung, 1987). 


A brief example is proved here to illustrate the ability of the SI method to correctly model a 
complex, time variant system. The soil system being modeled is between 47 and 17 meters deep 
beneath the A array. The input signal is vdha47_16e, and the actual output is vdhal7_l6e. The 


system is being modeled as a 3-DOF system 


The first model used was a [6 6 1] ARX model. In this case we are making the assumption that the 
signals are stationary and a time-constant model can capture the mechanical behavior or the system 
being excited by event 16. Figure 4.1 shows the comparison between the actual output from the 
30 meters of soil compared to the simulated output. It is seen that the model does very well (RMR 
error = 1.11) for most of the duration of strong shaking, leading to the conclusion that we have 
captured the soil behavior with the stationary model. However, the model has trouble matching 


system behavior for the first 5 seconds of strong shaking. 


As discussed in Chapter 3, the natural frequency, damping factor (% of critical) and participation 
factor can be calculated from the ARMA system identification. The respective values for these 


mechanical properties based on the stationary ARX estimation is given on Table 4.1. 


Table 4.1 System properties based on ARX identification. 


ed meee 
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Comparison of actual (solid) and ARX simulated (dashed) output at 17m. 
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Fig. 4.1 Comparison of output simulated by stationary algorithm to actual system output. 


The same soil system was then modeled using a recursive Kalman filter algorithm to see if the 
time-variable nature of the soil system could be captured. The results of the RARX simulation is 
given in Fig. 4.2. It is seen that the recursive simulation captured all aspects of the time history of 
the system very well (RMS error = 0.58). The amazing congruity of the simulated and actual 
system output attests to the fact that we can capture the mechanical behavior of the system with a 


few (6) simple parameters. 


Comparison of actual (solid) and RARX simulated (dashed) output at 17m. 
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Fig. 4.2 Comparison of output simulated by recursive (nonstationary) algorithm to actual system 
output. 
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CHAPTER 5 - CONCLUSIONS 


For the initial interpretation work it was decided to concentrate on the events for which accurate 
pore pressure records exist and show some increase in pore pressure during strong shaking. Events 
12, 16, and 17 meet this criteria. A summary of the properties of these three temblors is given in 
Table 2. Events 12 (30 July, 1986) and 16 (14 Nov., 1986) were major events and have been 
discussed in detail (e.g. EPRI, 1989; Chang et al., 1991; Anderson, 1993). Event 17 is an 
aftershock of event 16 which generated enough pore pressure to register on the recording 
equipment. The time histories for acceleration, velocity, and displacement are shown in Fig. 3. 
The corresponding traces for the largest acceleration components for events 16 and 17 are shown 


in Fig. 4 and 5. 


Event 7 is considered to be similar to event 12 (Elgamal, 1994), however the pore pressure 
transducers were not installed at this time. Event 7 is also important since it was analysed by 
previous researchers (Chang et al., 1989, 1990, 1991), so further analysis will build on this base of 
work. We will make use of this fine work by evaluating the statistical similarity of events 7 and 
12. If there is stochastic congruency between the two, the evaluation of event 12 pore pressure data 


will further our understanding of the site response to strong motion shaking. 


We are now ready to continue the application of SI to the Lotung data. If previous results are any 
indication, a great deal about the nonlinear behavior of earthquake excited soils will be learned 


(Glaser, 1995a). 
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APPENDIX A: COMPLETE SET OF VELOCITY TIME HISTORIES FOR EVENT 16, 
LOCATION A. 
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APPENDIX B: PORE WATER PRESURE TIME HISTORIES FOR EVENTS 12, 16, 17. 
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APPENDIX C - MATLAB PROCEDURE TO INTEGRATE ACCELERATION RECORDS TO 
VELOCITY AND DISPLACEMENT. 


function [a,v,d]=vd(f,Lcut,dT,n,pre,name) 

% function [a,v,d]=vd(f,Lcut,dT,n,pre) 

% 

% detrends and filters the input acceleration, 

% and integrates twice to give velocity and displacement. 

% The results are plotted so they can be reviewed. 

% 

% f is the input acceleration vector 

% Leut is the low-end cutoff frequency in Hz. 

% AT is the time step 

% n is the order of the Butterworth filter; 

% pre is the pre-event segment length to be zeroed. 

% 

% since the acceleration is filtered twice, the 

% effective order of the filter is double the value of n. 

% 

zip=(1:pre); 

zip=zeros(size(zip)); 

% 

[b,c]=butter(n, Lcut*dT*2.0, ’high’); 

a=filtfilt(b,c,f); 

a=dtrend(a(:,1),1,pre); 

a=detrend(a); 

a=a-a(pre+1); 

a(1:pre)=zip; 

% 

v=inttrap(a,dT); 

v=filtfilt(b,c,v); 

v=dtrend(v(1,:)’,1,pre); 

v=detrend(v); 

v=v-v(pre+1); 

v(1:pre)=zip; 

% 

d=inttrap(v,dT); 

[b,c]=butter(n, (Lcut*dT*2)/1.0, *high’); 
=filtfilt(b,c,d)’; d=dtrend(d(:,1),1,pre); 

d=detrend(d); 

d=d-d(pre+1); 

d(1:pre)=zip; 

% 

time=(length(f)*dT)/1.0; 

t=[dT: dT: time]; 

temp=name(1:(length(name)-6)); 

clf 

subplot(3,1,1), plot(t,a(1:length(t))), title([temp,’: acceleration’]), grid on 

subplot(3,1,2), plot(t,v(1:length(t))), title({temp,’: velocity’]), grid on 

subplot(3,1,3), plot(t,d(1:length(t))), title({temp,’: displacement’]), grid on, xlabel(’Time (seconds)’) 
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